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Introduction 


There is a need to derive a power spectral density (PSD) envelope for nonstationary acceleration time 
histories, including launch vehicle data, so that components can be designed and tested accordingly. 

Three methods are considered in this paper using an actual flight accelerometer record. 

The first method divides the accelerometer data into segments which are idealized as "piecewise 
stationary" in terms of their respective PSDs. A maximum envelope is then drawn for the superposition 
of segment PSDs. This method initially requires no assumptions about the response characteristics of the 
test item, but vibration response spectra may used for peak clipping as shown in the example. 

The following two methods apply the time history as a base input to a single -degree-of-freedom system 
with variable natural frequency and amplification factors. The response of each system is then calculated. 
Upper and lower estimates of the amplification factor can be used to cover uncertainty. 

The first of this pair is the energy response spectrum (ERS), which gives energy/mass vs. natural 
frequency, as calculated from the relative response parameters. 

The final method is the fatigue damage spectrum (EDS), which gives a Miners-type relative fatigue 
damage index vs. natural frequency based on the response and an assumed fatigue exponent, or upper and 
lower estimates of the exponent. 

The enveloping for each of the response spectra methods is then justified using a comparison of candidate 
PSD spectra with the measured time history spectra. The PSD envelope can be optimized by choosing the 
one with the least overall level which still envelops the accelerometer data spectra, or which minimizes 
the response spectra error. 

This paper presents the results of the three methods for an actual flight accelerometer record. Guidelines 
are given for the application of each method to nonstationary data. The method can be extended to other 
scenarios, including transportation vibration. 
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Sample Flight Data 


FLIGHT ACCELEROMETER DATA - SUBORBITAL LAUNCH VEHICLE 
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Figure 1. 


Derive an optimized PSD with 60-second duration to envelope the nonstationary, base input 
acceleration time history in Figure 1, via each of the three methods. This PSD is needed for a 
hypothetical avionics component which is to be mounted near the measurement location on a 
future flight. The component must be designed and tested accordingly. 

An option is to cover the liftoff event from 0 to 2 seconds as a separate shock or transient 
vibration event, but the entire record will be used as a test of the methods. 

Furthermore, a ramp-plateau- ramp envelope PSD will be derived in each of the three cases. 


Piecewise Stationary Method 

There are numerous ways to implement the piecewise method. 

For this example, the data is divided into consecutive 2.5-second segments. Each segment is idealized as 
being stationary. A narrowband PSD is then calculated for each segment. A maximum envelope PSD is 
taken for the superposition of segment PSDs. The envelope is converted to 1/12 octave format, as shown 
in Figure 2. 

The 1/12 octave curve is then enveloped by the Simplified PSD as shown in Figure 2. Simplification is 
advisable because the frequencies of the spectral peaks may potentially have flight -to- flight variation. The 
enveloping is performed using vibration response spectra per Reference 1. This method allows for peak 
clipping. A single amplification factor of Q=10 was used for the enveloping, but a pair of Q values could 
be used for added rigor. 
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POWER SPECTRAL DENSITY - PIECEWISE ENVELOPE CASE 
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Figure 2. 


Energy Response Spectrum Method 

The ERS method is taken from Reference 2. The concept seems to have originated in earthquake 
engineering and has since been a topic of research by Sandia National Labs. The equations for 
this method are shown in Appendices A & B. There are three energy components: kinetic, 
dissipated, and absorbed. The energy terms are calculated from the relative response terms of 
the spring- mass system with variable natural frequency and amplification factor. 

The ERS could hypothetically be applied by taking any of these components or some 
combination thereof. The approach in this example is to take the summation of each of the three, 
which yields the total input energy. Two amplification factors are used with Q=10 & 30. 
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POWER SPECTRAL DENSITY - ERS FLIGHT DATA ENVELOPE 2.3 GRMS 
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Figure 3. 


The method was applied by generated a series of candidate PSD functions. A corresponding 
time history was then synthesized for each PSD. 1 The ERS was calculated for each time history. 
The PSD and its ERS were then scaled so that each ERS would be greater than or equal to that of 
the flight time history at each natural frequency, with the two ERS curve equal at least one 
frequency. The winning PSD was that which yielded the least overall ERS error for the 
combined Q cases. 

The resulting optimized PSD for the flight data is shown in Figure 3. The ERS justification for 
this envelope PSD is shown for the two amplification factor cases in Figures 4 & 5. The Q=30 
case drives the envelope PSD. 


1 A method for directly calculating an ERS from a PSD is a topic for a future paper. 
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ENERGY RESPONSE SPECTRA Q=10 
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Figure 4. 


ENERGY RESPONSE SPECTRA Q=30 
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Figure 5. 
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Fatigue Damage Spectrum Method 


The response of the spring-mass oscillators to the flight accelerometer time history can be calculated in 
the time domain via the method in Reference 3. The rainflow cycles can then be calculated using the 
method in Reference 4. A relative damage index is calculated using a Miners -type summation for each 
permutation of natural frequency, amplification factor and fatigue exponent. 

The process can be carried out in the frequency domain from the response PSD for each candidate input 
PSD. The Dirlik method is then used for the FDS calculation, as taken from Reference 5 & 6. 

The optimization for this case was done by selecting the candidate PSD which had the least overall levels 
for displacement, velocity and acceleration. Another possibility would have been to select the one which 
gave the least error in terms of the FDS comparisons. Both approaches would mostly likely yield similar 
envelope PSDs, but verification is pending. 

Further details for these calculations are given in Appendices C & D, and in Reference 7. 

The amplification factors for the example are Q=10 & Q=30. The fatigue exponents are b=4.0 & 6.4. 
This gives four combinations. The selected fatigue exponents are those common for avionics 
components. 

The resulting optimized PSD is shown in Figure 6. 


POWER SPECTRAL DENSITY - FDS FLIGHT DATA ENVELOPE 2.4 GRMS 
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Figure 6. 
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Legend: Flight Data 


PSD Envelope 


FATIGUE DAMAGE SPECTRA Q=10 b=4 


FATIGUE DAMAGE SPECTRA Q=10 b=6.4 
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FATIGUE DAMAGE SPECTRA Q=30 b=4 
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FATIGUE DAMAGE SPECTRA Q=30 b=6.4 
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Figure 7. 


The corresponding FDS plots are shown in Figure 7. The case with Q=30 and b=6.4 is the driver for the 
envelope PSD. 
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Envelope Comparison 
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Figure 8. 


The PSD envelope comparison is shown in Figure 8. 

The Piecewise PSD is the most conservative in terms of overall level, as expected. The FDS -derived PSD 
shows the justifiable reduction within the scope of the assumed amplification factor and fatigue exponent 
combinations. 

The ERS-derived PSD is the least of the three in terms of overall level. 


Conclusion 


Selecting the enveloping method for a given nonstationary time history is a matter of engineering 
judgment. 

The piecewise method is the traditional approach and requires the fewest assumptions, but it may be 
overly conservative for many cases. 
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The FDS method requires the most assumptions , but appears to be most sound because fatigue is the most 
likely failure mode for random vibration environments. 

The ERS method needs justification that accumulated energy/mass correlates with damage in aerospace 
components. 2 

The FDS method seems the most sensible of the three, but further research is needed. 


Post Script 

The Matlab scripts for performing the calculations in this paper are given at: 

https y/vibrationdata.wordpress.com/2013/05/29/vibrationdata-matlab-signal-ana lysis -package/ 
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2 Note that peak energy is typically used in earthquake engineering in a roundabout way such 
that the resulting spectrum is really a pseudo velocity spectrum 
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APPENDIX A 


Eg uatio n o f Mo tio n 

Consider a single degree-of- freedom system subjected to base excitation. 


k 



y(t) 


Figure A-l. 


where 


m 

is the mass 

c 

is the viscous damping coefficient 

k 

is the stiffness 

X 

is the absolute displacement of the mass 

y 

is the base input displacement 


The double-dot denotes acceleration. 


The free-body diagram for the system is 
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Figure A- 2. 
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Summation of forces in the vertical direction, 


^F= mx (A-l) 

mx= f (t) + c(y — x) + k(y - x) (A-2) 

Define a relative displacement 

z=x-y (A-3) 

Substituting the relative displacement terms into equation (A-2) yields 

m(z + y) = — cz — kz (A- 4) 

mz + cz + kz = -my (A-5) 

Dividing through by mass yields, 

z + (c/m)z + (k/m)z = — y (A-6) 

By convention, 

(c/m) - 2^co n (A-7) 

(k/m) = co n ^ (A- 8) 

where co n is the natural frequency in (radians/sec), and c is the damping ratio. 

Substituting the convention terms into equation (A-6) yields 

9 

z + 2£,co n z + co n “z = — y (A- 9) 

The equation of motion can solved by the method in Reference 3 for the case of an arbitrary base 
input. 

Furthermore, allow the natural frequency to be an independent variable. The task is to calculate 
the “damage” at each natural frequency of interest for a fixed Q value. 

Q = 1/(24) (A- 10) 
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APPENDIX B 


Energy Equations 

The input energy per mass Ej is defined by integrating the equation of motion over the relative 
displacement. 


Ej = J q Z z(t)dz + 2£,co n J Z z(t)dz + co n 2 J Z z(t)dz 


(B- 1) 


The kinetic energy per mass Eg is 


E k = j*z(t)dz = £ z(t)z(t)dt 


(B-2) 


The dissipated energy per mass Ep> is 


E d = 2£,co n J Z z(t)dz - 2£,co n z(t) 2 dt 


(B-3) 


The absorbed energy per mass Ea is 


e a = co n 2 j z z(t)dz = m n 2 £ mm dt 


(B-4) 


Note that the energy terms are functions of the relative response amplitudes. 
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APPENDIX C 


A relative damage index D can be calculated using 


D = 



i=l 


(C-l) 


where 

A- is the response amplitude from the rainflow analysis 
n. is the corresponding number of cycles 
b is the fatigue exponent 


Note that the amplitude convention for this paper is (peak- valley)/2. 
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APPENDIX D 


Dirlik Method 


The nth spectral moment m n for a PSD is 

m n = J o °°f n G(f)df 


(D-l) 


Where 

f is frequency 

G(f) is the one-sided PSD 


The expected peak rate E[P] is 


E[P] = V m 4 /m 2 


(P-2) 


The Dirlik histogram formula N(S) for stress cycles ranges is 


N(S) = E[P J • T • p (S) 


(D-3) 


where 


T is the duration 

S is the stress cycle range (peak-to-peak) 


The function p(S) is 


P (S) = 


Dp 
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(D-4) 
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The coefficients and variables are 


1 + y 2 


l-y-D^ +D^ 

Do = 

2 1-R 

(D-6) 

D 3 =1-D! -D 2 

(D-7) 

z- s 

(D-8) 

1.25(y-D 3 -D 2 R) 

Dl 

(D-9) 

y-Xm-D 2 

— o 

l-y-Di + 

(D-10) 

m 2 

V m 0 m 4 

(D-ll) 

mi m 2 
x m — ■> 

m 0 V m 4 

(D-12) 


A cumulative histogram of the peaks can then be calculated from equation (D-3). 

The stress range of individual cycles can then be obtained by interpolating the cumulative 
histogram 
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